#ifndef __MESH__
#define __MESH__

#define EPSDIST (1.0e-4)


//
// 1. data structure:
// linked sorted list are used for storing all face infomation
//
typedef struct face{
  int _id;
  int _vtx[3];   // three vertex
  int _bnd[3];   // three edges
  int _nei[2];   // left and right tet
  int _boundary_mark;  // used as a mark for the boundary
  struct face *_next;
} Face;


// in Original mesh format
void readMeshHeadA(char *meshPath,
		  int *N_node, int *N_edge, int *N_element,
		  int *N_point, int *N_line, int *N_triangle);



void readMeshDataA(char *meshPath,
		  double *Node,    int N_node,     //
		  int    *Edge,    int N_edge,     //
		  int    *Element, int N_element,  //
		  int    *Point,   int N_point,    //
		  int    *Line,    int N_line,
		  int    *Triangle,int N_triangle,
		  int    *Triangle_Element,   // for find neighbour
		  int    *Element_Edge,        // edge-element 
		  int    *Triangle_Edge,       // boundary triangle
		  int    *Line_Edge,
		  int    verbose);           // boundary dof

// write in Gmsh format
void writeGmsh(char *meshFileName,
	       double *Node,     int N_node,
	       int    *Element,  
	       int    *ElementMark, int N_element,  
	       Face   TriangleList, int N_triangle,   
	       int    verbos);

// read mesh in current format
void readMeshHead(char *meshFileName,
		  int  *N_node,
		  int  *N_edge,
		  int  *N_triangle,
		  int  *N_element,
		  int   verbose);

void readMeshData(char *meshFileName,
		  double *Node,    int N_node, 
		  int    *Edge,    int N_edge,
		  int    *Element,        // node dof for the element
		  int    *ElementEdge,    // edge dof for the element
		  int    *ElementMark,    // ElementMark for the element
		  int     N_element,  
		  Face   *TriangleList, int N_triangle,
		  int verbos);           // boundary dof

// write mesh in new format
void writeMesh(char   *meshFileName,
	       double *Node,    int N_node, 
	       int    *Edge,    int N_edge,
	       int    *Element,        // node dof for the element
	       int    *ElementEdge,    // edge dof for the element
	       int    *ElementMark,    // ElementMark for the element
	       int     N_element,  
	       Face   TriangleList, int N_triangle,   
	       int verbos);           // boundary dof

// read from AFEPack mesh
void readAFEPackMesh(char *meshFileName,
		     double **Node,    int *N_node,     //
		     int    **Element, int *N_element,
		     int verbose);

// tool function to build mesh boundary and edge form incomplete mesh data
// Warning: THIS IS TIME CONSUMING!
void build_mesh_boundary(int **Triangle, int **TriangleEdge, 
			 int **TriangleElement, int *N_triangle,
			 int *Element, int N_element,
			 int *Edge,    int N_edge,
			 int verbose);

void build_element_edge(int *ElementEdge,  // output
			int *Element, int N_element,
			int *Edge, int N_edge,
			int verbose);

/**
 * read the gmsh file only with nodes' and tetrahedrons' informaiton
 * including flags(boundary mark)
 */

int readGmshHead(char *meshFileName,
		 int  *N_node,
		 int  *N_triangle, int *n_flag_per_tri,
		 int  *N_element, int *n_flag_per_elem,
		 int   verbose);

int readGmshData(char *meshFileName,
		 double  *Node, int N_node,
		 int *Triangle, int N_triangle,
		 int *TriangleFlag, int n_flag_per_tri,
		 int  *Element, int N_element,
		 int *ElementFlag,  int n_flag_per_elem,
		 int   verbose);

/**
 * write out mesh data according to the ElementFlag
 * only Nodes' position and Element vertex are output in Gmsh Format
 */

void write_gmsh_marked(char *mesh_prefix, 
		      double *Node, int N_node,
		      int *Element, int N_element,
		      int *ElementFlag, int n_flag_per_elem, 
		       int flag_idx);

/**
 * If the index of any tetrahedron is not in positive order
 * change the first and the second nodes, which yielding a positive order
 * The output of Gmsh should obey this rule and nothing has to be changed
 */
double volume_tet(double *v0, double *v1, double *v2, double *v3);

void reorder_node(double *Node, int *Element, int N_element);

void build_edge(int *Edge, int *N_edge, int *ElementEdge,   // output
		int *Element, int N_element, int N_node);     // input mesh

Face * find_face(Face *FaceList, int i, int j, int k);

void build_face(Face *FaceList,  int *N_face,   //output, insert all faces
		int *Element, int *ElementEdge, int N_element, // input
		int *Edge,    int N_edge, int N_node);

void build_triangle(Face *TriangleList, int *N_triangle, //output
		    Face *FaceList, int *ElementMark, int N_node);

double element_volume(double *v0, double *v1, double *v2, double *v3);
double volume(double *Node, int *Element, int N_element);

void mirror_mesh(double **Node, int *N_node, 
		 int **Element, int *N_elem, int axis);

int 
load_mesh_info(char *fileName, int *N_node, int *N_element);

int 
load_mesh_data(char *fileName, 
	       double *Node, int N_node, 
	       int *Element, int N_element);

// write a real-valued field into DX format
void writeMeshDX(char *fileName,
		 double *Node, int N_node,
		 int *Element, int N_element);  


#endif
